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ABSTRACT 

There are a variety of meteorological forecast problems which require 
high spatial resolution in only a limited area* An important example of 
this type of problem is the prediction of tropical cyclones. This study 
tests a simple finite element prediction model with a variable element 
size. The shallow water equations are used and the motion is confined in 
a periodic channel on a f-plane. The Galerkin technique is applied to 
linear basis functions on triangular elements. The model uses leapfrog 
time differencing and periodic restarts. The model is tested with a wave 
imbedded in a mean flow and also with an isolated vortex. The experiments 
with a uniform element size show excellent phase propagation, but some 
small scale noise is generated. The introduction of momentum diffusion 
terms helps to control the noise. The model is also tested with elements 
which decrease abruptly in scale along a line and with elements which de¬ 
crease smoothly. Both of these cases generate more noise than with 
uniform elements. 
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I. INTRODUCTION 


Operational nijmerical models, because of computational limitations, 
are able to depict and predict atmospheric processes only down to a 
limited scale. Until the introduction of the National Meteorological 
Center's limited fine mesh (LFM) model, which is superimposed on the 
coarse mesh Northern Hemisphere stereographic grid, grids were uniform 
on any given map projection and the advent of LFM has brought increased 
accuracy in wave phase and amplitude prediction due to smaller trunca¬ 
tion errors [Houghton and Irvine (1976)]. Nevertheless, problems arise 
at the boundaries of the fine mesh grid because of the abrupt change in 
grid size and the need to furnish boundary values interpolated from the 
coarse tnesh. Hence it seems logical that a scheme which permits a 
gradual change in resolution would be desirable. 

A relatively new technique, known as the finite element method (FEM), 
seemed suited to the problem of a non-uniform grid. In theory the de¬ 
pendent variable may be predicted by this method at any number of loca¬ 
tions chosen at will. The main thrust of this particular investigation 
was therefore to test this hypothesis. Secondary emphasis was placed 
on determining the computational limitations, computer time and storage 
requirements of such a procedure. 

Only recently have the atmospheric sciences seen the use of the 
technique. Among the notable contributions are those of Cullen (1973, 
1974, 1976); Lee, Gresho, and Sani (personal communication); and Hinsman 
(1975). Since the finite element method (FEM) is an implicit scheme in 
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space, physically reasonable solutions involve solving a large system 
of linear equations. This implies substantial memory requirements and, 
depending on the technique of matrix inversion, possibly considerable 
CPU time. 

Specifically, the technique of Galerkin (1915) was used to transform 
the equation set into an integrable form. Galerkin*s method has recent¬ 
ly been proven invaluable in engineering disciplines. For example, the 
field of indeterminate structural analysis, pressure vessel and air 
frame design, and continuum mechanics have all greatly benefited. 

The equation set used was the "shallow water”, barotropic equations 
which approximately describe atmospheric motion. These equations model a 
simple, physically reasonable form of atmospheric motion. Additional 
simplification was gained by assuming the coriolous parameter, f, was 
constant and equal to its value at the mid-channel latitude. This assump¬ 
tion does not degrade model performance significantly since the actual 
value of f varies slowly over the domain studied. The value of this 
approximation lies principally in reducing the number of terms in the 
equations by two and eliminating the "Beta effect”. With no f variation, 
sinusoidal waves can move in an east-west direction without any resultant 
latitudinal tilt. 

Two basic types of domains and initial conditions were tested. Both 
domains were square with a channel length of approximately 3500 kilometers 
The channel width was equivalent to a region extending from the equator 
to 30 degrees north latitude. The channel was periodic in the east-west 
direction. 
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Simple experiments were carried out for sinusoidal initial conditions 
which consisted of a single wave in the longitudinal or x-direction and 
a half wave in the latitudinal or y-direction, A geopotential surface 
tilt was added in the y-direction to allow for any mean zonal flow imposed. 
Numerous tests were performed with these initial conditions using various 
domain subdivisions. All domain subdivisions were done in a regular manner 
by varying the grid spacing in either or both the x and y directions over 
the domain. Additional tests were also conducted on some of these meshes 
with an analytic vortex. 

Final and more’ complex tests involved a domain on which the grid reso¬ 
lution decreased regularly toward the center of the domain. Initial condi¬ 
tions used consisted of those mentioned above and a localized disturbance 
or vortex which spanned the central region of varying resolution. 

The simpler tests were undertaken to determine the effects of various 
smooth and abrupt grid changes on model performance. The final, complex 
model was evaluated with an eye towards ultimate implementation as an 
operational tool for forecasting tropical cyclone movement. It was hoped 
that the flexibility of choosing a varying grid resolution would be a 
significant improvement over current finite difference models. It seems 
desirable to have a scheme by which one can vary grid resolution smoothly 
to accommodate the scale of motion changes that occur proceeding radially 
inward into a typhoon. 
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II. FINITE LINEAR ELEMENTS 


As is true with the spectral methcxi, the FEM assumes that the depen¬ 
dent variable can be represented by a known function of space. The two 
techniques differ since spectral approximations are continuous and are 
infinitely differentiable while FEM functions are defined very locally 
(over an element)^ vanishing elsewhere, with limited continuity. Gener¬ 
ally, FEM functions are low order polynomials, defined as planar surfaces 
locally in the case of linear elements. That is, for a given variable, 

(p , one assiomes that over each element the basis or local support function 
has the form: = a + bx + cy where a, b, and c are constants. 

Figure 1 shows an exploded, three dimensional view of a small domain 
which has been subdivided into triangular elements. Above this domain 
lies the group of linear basis functions which are pyramids rising to an 
approximation of (p over a nodal point. Above these pyramids is the seg¬ 
mented, planar approximation of (p with dashed cross-sections of the exact 
(p surface at evenly spaced Ax values. Note that <p is shown varying 
only in one direction in this case and that the (p approximations, $ , 
are shown to be exact at the nodal points. The latter fact is, of 
course, not true in general. What is required for piecewise continuity 
in <p^ is that no nodal point may be located on the side of a triangle. 

Once a form for the unknown as been assumed, Galerkin*s technique may 
be applied. In its basic form, this method assumes an unknown has the 
following approximate form locally over a nodal point, j and its surround¬ 
ing elements; 

= YjV^ (II-l) 
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Figure 1. Linear finite element approximation. 
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(with the repeated svibscript implying a sum over that subscript) 


Vj is defined as the basis or support function and has maximum value 
of unity at the node^ sloping to zero at surrounding elements* outer 
boundaries. Thus with the aid of Figure 1 it is seen the Yj = 4^j 
any nodal point j• 

To solve a differential equation in (}> , assume (p is of the form 
(II-l), multiply by a suitable "test” function, and integrate over the 
domain. These test functions, if chosen wisely, greatly simplify the 
resulting equation and make it easily integrable. A generalized example 
treats, in one dimension, the equation (with ^ { } a differential 

operator): 


c^{(p} - f (x) =0 0 ^ X ^ a 


An approximate form for (p is assxamed, (p ; the test function, 
g(x), is applied, and the integration performed. The resulting equation 
is: 

a 

J* [a!<^{(p} - f(x)]g(x) dx = 0 (II-2) 

o 


Since [ <^{(j)}-f (x) ] = R where R is the error in the approximation, 
this integration should minimize R, with g(x) properly chosen. Galerkin's 
choice for g(x) is the same form as that for (p^ , (II-l) , with the co¬ 
efficient Yj E 1. The attraction of this approach is that since both 


g^ and (p^ are defined only locally, i.e. over elements surrounding 
nodal point j, the result is a linear matrix equation which is banded. 

This results from a partial integration of (II-2). Another useful feature 
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is that linear elements have domain integrations which are analytically 
obtainable. 

This direct integration is made straightforward through the use of 
a natural or area coordinate scheme for each element as explained in 
Desai and Abel (1972)• In two space dimensions the location of a point 
P(Xry) inside an element (see Figiire 2) can be defined by describing the 
three areas: A^^, A^/ and A^, subtended by lines from P to the triangle's 
vertices (x^, y^, j=lf 2, 3), Zienkiewicz (1971) shows how the location 
of P can be transforraed into area coordinates via: 


X = L.x. 

3 3 

y = L.y. 

1 - Li + t L3 


(II-3) 


with; 


j = 1, 2, 3 


L.= A./A 
3 3 

A = total elementary area 


If the differential equation to be solved is; 


94) 3(t) , . . 

*1 i i *2 


= 0 


(II-4) 


1^2 


Galerkin's technique yields: 


■'2 2 


Y 


. / / (V. + V. + V.)V.dxdy 

2 . J J IX ly 1 j 


= 0 


^1 ^1 


with variable subscript implying partial differentiation. 
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Figure 2. Area coordinates. 
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This matrix equation may be solved for the coefficients once 

the area integration is .performed. Zienkiewicz show that for triangular 
elements with linear basis and test functions a general area integration 
formula is: 



_ m _ n , - 

dxdy = 


mini 
(m+n+2)! 


2A 


(II-5) 


It is noted that, for m, n>l, non-linear terms are being formulated. 

A 

For example, with m=l and n=2, the value of (II-5) becomes — . 
Transforming equation set (II-3) into matrix form and solving for 
tlie L vector one obtains: 



with a^, b^ as defined in Figure 3. 


Differentiation of (II-6) shows the differential operations desired 


to complete the integration of (4) are: 


3x 2A 3 l. 

1 

A. = -LA. 

3y 2A 3 l. 

1 


(Henceforth a repeated subscript will imply a sum from one through 
three.) 
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Figure 3. Local coordinate system. 
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Applying the chain rule in natural coordinates: 


3 V. 


3 V. b^ 

3 V. b, 3 v. 

1 

_1 

1 , 2 

—L + ^ 1 

9x 

2A 

2 A 

3L2 2 A 3L2 


but 


3V. 

1 



and since 


0 for i ?£ k 
3V. 

V. = L. ; - 5 ^ =1 for i = k 

3. X oXj- 


3V. 

b. 


3V. 

a. 

1 

_ 1 

and 

_ 1 

_ 1 

3x 

2A 

3y 

2A 


(II-7) 


Therefore equation (II-4) is easily integrated on an elementary level. 
These results can then be distributed throughout the resultant coeffi- 

/S 

cient matrix and the solution vector {({)} obtained via a linear equation 
solver. 
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III. EQUATIONS AND BOUNDARY CONDITIONS 


A form of the barotropic, shallow water, primitive equations was 
selected which was developed by Phillips (1959). The original intent 
of this study was to develop a complete, barotropic model for a tropical 
channel. Time limitations forced the use of a model which was purely 
cartesian. That is, although the original equation set was designed 
in mercator coordinates, a constant map factor of unity was used. In 
theory, and based on other work done by this author, extension of the 
Galerkin formulation of the equation set to include complete sphericity 
is easily done. 


The set in cartesian coordinates becomes 

11=<''♦'1 




3u ^ 8 u, . . Q 

+ u + V + 2fiv sin 0 
9t 9x 9x 9y o 


■ 5 — = - + u -^ + V •?—] - 2fiu sin 0 

ot oy Ox 9y o 


(III-l) 

(III-2) 

(III-3) 


with the transform from spherical to the cartesian system accomplished 
using; 


= aX 


and 


= a In 


cos 0 
1 -sin 0 

latitude 

longitude 


0 = 15“ 

o 
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Ihe only mathematical boundary conditions imposed were those of no cross 

channel flow at the latitudinal boundaries and of cyclic continuity in 

the longitudinal direction. Through the course of model development 

other conditions were tried at the boundaries due to physical constraints. 

A geostrophic balance at channel walls was imposed in the mass continuity 
1 

equation(III-l ), i.e.: 


3y 


2Uu sin 0 

o 


(III-4) 


Additionally it was thought necessary to impose a zero normal 
gradient of zonal wind at the walls. During some early experiments 
with an abruptly varying element size, it was found that imposing this 
condition on the zonal equation of motion yielded instability after six 
to 12 hours of time integration. Later experimentation showed the zonal 
equation required no explicit boundary conditions. 


This condition is consistent with that of no cross boundary flow 
and follows from inserting V=0 in equation (III-3). 
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IV. DIFFUSION 


Early experiments with a coarse mesh indicated that some type of 
filtering would be required since short waves of twice the coarse grid 
size were readily generated. These short waves were evident in the dis¬ 
turbance patterns of the geopotential and zonal wind fields. Since ex¬ 
perience and knowledge of atmospheric FEM applications is still in its 
infancy it is not known if this phenomenon is an artifact of the time 
stepping scheme employed, the boundary conditions imposed, or the basic 
technique. 

A simplistic, second order diffusive damping term was included in 
the equations of motion. To keep the Galerkin formulation of the 
diffusion consistent with first order terms, an area integration was 
performed and Gauss' divergence theorem was applied. Consider the diffu¬ 
sion term in the zonal equation (III-2) with Galerkin*s technique applied 


K y'V^uV.dA = K f V'VuV.dA 


H 


H 


K„ {>• V. Vu*n dr -S VV. •VudA} 
“ A ^ 


(IV-1) 


where dr is differential distance along the path of integration on the 

/V 

outer boundary, K is the horizontal viscosity coefficient and n is 
H 

the unit normal to the area or domain, A. 


The contour integral in (IV-1) vanishes due to cyclic continuity and 
9u 

an implied = 0 . The diffusion term in the meridional equation is 
advanced similarly. 
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V. EQUATION FORMATION 


As an example of the use of the Galerkin technique consider equation 
(III-l) . A consistent form for (J) was assumed and substituted into 
(III-l). Then an area integration was performed. 

Defining the inner product of a basis function g(x,y) and a test 
f\anction, V(x,y) as: 

/ / gV dxdy = <g,V> 

y X 

An application of the FEM to the continuity equation gives: 

"■!? <“*>- " w <vf), v> 

A partial integration of the x-derivative terra yields: 





V> 



(u(t)V)dA - 


/u*|^dA 


But the first term on the right hand side vanishes because of cyclic 
continuity in x. Similarly, since v = 0 on the boundary, the correspond¬ 
ing term in the partial integration of the y-derivative term also 
vanishes. 

Applying Galerkin's linear technique it is assumed (j), u, and v are 
of the form 


4 ) = y.V. 

D 1 


u = a.V. 
D 3 

V = 6^v. 
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which, with a concurrent test function assvunption, transforms the con¬ 


tinuity equation to: 


<V. V.> = 
1 


Y. [<a, V, V. V. > 
j k k 3, IX 




where the dot denotes a time derivative.’ 

Using a centered finite time difference the final form of the 
equation is: 


v.> = 2At yn [<a, ‘V, V. V. > + <$”v, V. v. >] (v-l) 

3 3 3f 1 '] k k ], IX k k 3. IX 


with h being a time level. 

Equation (V-l) is a matrix equation of the form 


[A] {x} = {b} 


where; 


[A] = <V. V.> 

3' 1 




and 


{b} the right-hand side of (V-l). 


23 



1 : : 






The u and v momentum equations were similarly advanced including 


diffusion and become, respectively: 


-a.’^[<a.Vv. ,v.>+<B "v, V. ,v.>] 
3 k k jx' 1 ^k k dy' 1 ' 


(a 




,v. >+<a.\ , V. >] 
H k kx IX k ky ly 




+ 2!)<sine^a,^\ , V.>+<YkV^y,V^> )(V-3) 

-■S, V 


At this point it is noted that all the inner products are functions 
only of space and need be computed but once. Also useful is the fact 
that the coefficient matrix, identical in all three equa¬ 

tions, Finally, at this stage all equations are uncoupled at any time 
step. 

As is apparent, the time discretization scheme chosen was the leap¬ 
frog or centered difference method. It was chosen primarily for its 
simplicity. Obvious drawbacks are its inherent computational mode and 
the somewhat more restrictive maximum time step than some semi-implicit 
schemes may allow. The latter factor was expected to and did play a 
major role in the amount of computational time involved since the 
smallest grid resolution essentially limits the maximum allowable time 
step. It was expected, based on the v;ork of Cullen (1973) using 
simple linear advection, that the time step criterion would be slightly 
more restrictive than a finite difference scheme with equivalent 
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m 


resolution. The leapfrog scheme used a half-forward, half leapfrog time 
step to start. Also, early experiments showed strong solution separation. 
This tendency was kept in check with a periodic restart, identical to 
the initial start, being applied every 12 time steps. 
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VI. TREATMENT OF BOUNDARY CONDITIONS 


The cyclic boundary condition is natiirally included in the system by 
defining how the nodal points are connected. Remaining conditions are 
those imposed at the northern and southern walls. The geostrophic condi 
tion is essentially independent of time, so, to fit this condition into 
the model, boundary equations at distinct times, n+1 and n-1# were sub¬ 
tracted. That is: 

-2nsin9 u) ^ -2Qsin0 u)^ ^ = 0 

dy o dy o 

The Galerkin formulation gives: 


(Y <v. ,v.> = -f ,v.> - <a.^”^v. ,v.>] 

D 3 DY 1 o j 1 : 1 


with 


f = 2nsin9 
o o 


(VI-1) 


The rigid wall boundaiy formulation merely required setting the 
Galerkin approximation for meridional wind to zero at the i boundary 
points initially, cind requiring no time change in these values: 

<v.,v > = 0 

J D J 1 

These boundary conditions were applied at all points, i, on northern 
and southern boundaries. 

It should be realized that the geostrophic boundary formulations 
imply that a different coefficient matrix is formed than that shown in 
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equations (V-1, V-2, V-2). The only terms in these matrices which are 
different are those in the boundary equations. 
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VII. COMPUTATIONAL TECHNIQUES 


A primary operational disadvantage of any finite element technique is . 
the requirement of large amounts of computer storage. This results from 
the system of linear equations. The problem is not as bad as it first 
seems since in any one equation of an N x N system a majority of coeffi¬ 
cient matrix elements are zero. The problem then reduces to inverting a 
system with a very sparse coefficient matrix. If one can utilize an 
equation solver which takes advantage of this property, core requirements 
can be minimized. 

Specifically, for the ith nodal point, the number of non-zero entries 
in the coefficient matrix for equation i is one plus the number of elements 
surrounding that nodal point. Some storage schemes take advantage of the 
banded nature of the coefficient matrix. The band width is determined 
by the largest different between the number, i, of the ith equation and 
the numbers of the nodal points locally connected to nodal point i. Mini¬ 
mal bandwidth is thereby gained by choosing an efficient nodal point 
numbering scheme. 

It was decided not to minimize bandwidth structure for any cases under 
consideration since the time involved in selecting an efficient number 
scheme for differing domains possessing cyclic continuity could be better 
spent. A disadvantage of a banded matrix scheme is the storage of some 
zero values; a possible problem when large systems are to be solved. This 
author's unbanded approach has the disadvantage of forcing the use of an 
iterative solution technique for the linear system. 
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To fully compact the coefficient matrix a novel technique was used 
which was developed by Salinas (private communication). Essentially, 
only non-zero coefficient matrix entries are retained. The coefficient 
matrix is stored as a linear vector. Three vectors herein called NUM, 
ISTART, and NAME, are required to use the method. NUM contains for 
each j nodal point the number of points interacting with j including 
itself. NAME gives the information of which nodal points interact with 
j including j. ISTART tells where, in NAME, the numbers of the points 
interacting with j begin. NUM and ISTART are both dimensioned N for a 
system with N unknowns. NAME is dimensioned to the length of non-zero 
entries in the coefficient matrix. That is, NAME'S dimension is the 
sum of the number of nodal point interactions. 

The matrix assemblage required is mathematically equivalent to 
multiplying each basis function on a .nodal point by nodal point basis 
by the correct test functions, and then integrating over the portion of 
the domain spanned by these functions. Graphically, as shown in Figure 
4, the basis function at j is shown outlined in heavy black over a 
group of elements around nodal point j. Portions of the test functions 
which interact with this basis function are also shown. 

Since, in the case of linear elements, the domain integration is 
analytically obtained over an element, the matrix assemblage of an inner 
product, such as <Vj,V^> , is done most efficiently on an elementary 
level. To accomplish this, a local numbering scheme is required for each 
element. This table was stored in the two dimensional array, NPTS (N,3). 
NPTS contains, for each element, e, a set of the three global nodal 
points which define e. For each element these numbers are stored 
sequentially according to number of the second dimension of NPTS. The 
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Figure 4* A basis function and portions of two test functions. 
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choice of the point, k, to start niimbering e was arbitrary, but, once 
chosen, the remaining two points were numbered in sequence proceeding 
in a positive sense around e. For example in Figure 4 element e is 
labeled NPTS (e,l) = k, NPTS (e,2) = 1, and NPTS (e,3) = j. 

To illustrate matrix assembly using an element by element technique 
consider again Figure 4. The computational technique required is for 
each point describing e; namely k, 1, and j; the inner product value be 
distributed to its proper location in the coefficient matrix, • 

Two of the three test functions over e are drawn labeled t^ and t^. 

Into the equation for point k a contribution from the interaction between 
b^ (not shown) and t^ must be added; that is, <b^,t^> (known analytically) 
must be added into the kth element of equation k. Also t^^'s interaction 
with b^ must be entered into the 1th element of equation k. Finally, 
tj's interaction with b^ must be added into the jth element of equation k. 
The same local dispensing of interactions must be carried out for the 
remaining nodal points, 1 and j of e» This procedure when done for all 
elements in the domain gives coefficient matrix of the 

equation. 

The inner products used may either be computed as needed during the 
equations * time integration or computed prior to time integration and 
stored for all time. The latter method was followed since the original 
computer system used was in IBM 360/67 which is relatively slow by today's 
standards but was able to handle the core requirements of the problems 
tested. Also, since varying size elements would be used and the equations 
require numerous time steps for even a six hour forecast, an inner product 
storage is most efficient. 
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VIII. DOMAINS CONSIDERED 


A. BASIC MESHES 

To gain experience with the FEM, numerous grid meshing schemes were 
considered. The first problem chosen was one which required a minimal 
amount of central memory storage and execution time. This mesh is pre¬ 
sented in Figure 5(a). Actual tests were conducted on grids similar to 
the remaining por^tions of Figure 5. In all cases the domain is a bounded 
channel, cyclic in the east-west direction with a width and height of 
3503 km. 

Higher resolution, regular grids were also used. One of these is 
shown in Figure 5(b). The domain was divided into 12 latitudinal and 
longitudinal subdivisions (12x12). A finer resolution used on this domain 
had twice the resolution as Figure 5(b), (24x24), and is not shown. A 

coarser scheme (not shown) had six latitudinal and longitudinal sub¬ 
divisions, (6x6). 

Irregular grids were also used. Figure 5(c) shows the first of these 
tested wherein a central region in the x-direction had half the external 
mesh size covering half the channel length. A similar grid was used with 
a varying resolution in y. These two were combined to give the grid 
shown in Figure 5(d). 

Other irregular grids were used in which, rather than abruptly chang¬ 
ing the resolution, a smooth gradation was employed. The x varying ver¬ 
sion is shown in Figure 5(e) along with the version which varies in x and 
y. Figure 5(f). This zonal and meridional variation was made with no 
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Figure 5. Regular and graded meshes 











































































































































particular plan other than a smooth transition to an inner region of 
half the external or near boundary element sizes. 

All the grids so far mentioned were coded to allow a reversal in 
diagonal slant. Intuitively/ from examination of Figure 5, it would 
be surprising if some bias was not introduced when domains are divided 
with a built in diagonal preference. These same grids were also tested 
with an alternating diagonal slant similar to that shown in Figure 6. 


B. FINAL MESH 

In order to study a relatively small-scale vortex, a meshing scheme 
was developed whose spatial resolution increased radially toward the 
center of the vortex. The grid developed is shown in Figure 6. The 
routine used to locate the nodal points was executed independently of 
the time integration and its results were read in prior to overall inte¬ 
gration. Input data included nodal point locations and the local 
elementary numbering scheme, NPTS. 

This program was made flexible in order to allow for a varying reso¬ 
lution over an even number of latitude/longitude blocks in the total 
domain. Also included was a scheme by which the radial resolution on 
the subdomain could be varied in a regular manner. This procedure involved 
decreasing the distance between nodal points along the diagonals which 
extend from the corner points of the subdomain by increasing even powers 
of an arbitrary factor, K. This is, for the five diagonal divisions shown 


of length d^: 


d. . (K) d„ 


with K < 1 d^ = ungraded domain diagonal length 

In the ungraded region the diagonal slants were varied as shown to minimize 

any bias. 
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Figure 6. Final mesh. 
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IX* INITIAL CONDITIONS 


A. SINUSOIDAL 

For the basic testing on the rectangularly subdivided grids such as 
shown in Figure 5, simple initial conditions were used which consisted of 
a single wave in the x direction and one half wave in the y direction* 

Such initial conditions which satisfy the non-linear balance equation with 
f constant are given below: 


d) = f A sina^ cosa« - f U (y-y ) + (b 
^ o 1 2 o m 

— 

u = U — cosa. cosa^ (IX-1) 

W 1 2 

27rA 

V = - —— sina^ sma 
L 

with A = arbitrary perturbation amplitude 

V7 = channel width = 3500 km 

X = channel length = 3500 km 
L 

U = mean zonal wind = 10 m sec ^ 
y^ = mid latitude value of y = W/2 

<b = mean free surface goepotential height chosen as 

^ 4 2-2 

4*9 X 10 m sec 

= TTyA^ 

= 2Trx/X^ 


Two values of A were tested; one which gave a maximum perturbation zonal 
wind of 5*5 msec ^ and another for 0*55 msec A mean zonal flow of 
10 msec ^ was imposed in some experiments with this initial condition. 
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The initial conditions with and without zonal flow are presented in 
Figures 7 and 8 respectively. The perturbation amplitude in both cases 
is for a perturbed zonal wind of 5.5 msec 


B. VORTEX 

Analytic initial conditions were used for a small scale perturbation. 
A tangential wind, was defined in cylindrical coordinates as: 


= A — sin^(cr) 
T r^ 


r < r 
— 0 


(IX-2) 


Vo, = 0 

T 


r > r 


with r = radial distance from center of perturbation 
r^ = maximiam radial extent of perturbation 
C = TT/r^^ 

A = arbitrarily chosen amplitude in msec ^ 

An analytic value of (J)*(r), the disturbance portion of geopotential, 
can be obtained from the gradient wind equation: 


V_ 


JL 

r 


o T 


9r 


(IX-3) 


Integrating (IX-3) radially outward from the center of the vortex 


(r=0) one obtains: 

4.' (r) = - ^[sin(2cr) - ^ sin(4cr) ] 


4c 


8 


-[cos(2cr) - cos(4cr)]} + 

1-1 ^ J-o 

8c 

fA r . , 1 X 1 

— -— sin(2cr)- :r cos(2cr)> 

4 4c - 2 

0 8c 

+ constant of integration 


(IX-4) 
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Initial sinusoidal conditions with U = 


10 msec 
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Figure 8. Initial sinusoidal conditions with U = 0. 
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The amplitude of the disturbance. A, was chosen to give a desired 
maximum tangential wind component of 25 meters per second, a realistic 
value for tropical models at 500 mb. The mean 500 mb geopotential 
height, <{)^, was added to this <J)* disturbance field. Also added was a 
mean slope field, similar to that in (IX-1), to allow for a mean zonal 
flow. 

The radius of the distiirbance, r^, was chosen so that the vortex 
spaneed the graded portion of the domain shown in Figure 6 equal to six 
largest grid distances; thus r^ was approximately 875 kilometers. The 
initial conditions are presented in Figure 9 with mean flow component 
shown only in the <j)’ field. 
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Figvire 9. Initial vortex conditions with U - 10 msec 
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X. RESULTS 


A. DISCUSSION 

1. Notation 

The use of various meshes on the same domain enabled a detailed 
study of the basic properties of the FEM. Obviously the order of case 
study proceeded from the coarsest, 6x6 grid mesh to an 18x18 mesh with 
varying resolution in both directions to the 24x24 mesh which had the 
best overall resolution. To simplify the task of referencing all cases 
considered (see Table I). 

All the above cases were integrated through 48 hours. A maximum 
stable time step was roughly found using Cullen's (1973) work as a guide. 
Cullen states that for simple, two dimensional, linear advection, a stable 
FEM time step is l/»/3 times that of the corresponding finite difference 
(leapfrog) scheme. For the 12x12 mesh with a mean flow of 10 meters per 
second, Cullen's criterion gives a stable time step of 8.8 minutes. A 
nine minute time step was found unstable on the 12x12 while an eight 
minute time step was stable. This verifies Cullen's analysis since a 
corresponding stable finite difference time step is slightly in excess 
of 15 minutes. Time truncation errors were not of primary concern in 
this study, so close to maximum time steps were always used. 

2. Harmonic Analysis 

A harmonic analysis routine of the geopotential fields was used 
at specified times in the integration. This analysis was accomplished 
along constant latitude circles for wave numbers through three, six and 
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TABLE I 


Experiment Abbreviations 
(all valid at 48 hours) 


6x6 

tSD] 

12x12 [SD] 

18x12 [SD] 

12x18 [SD] 

6x6 

[SID] 

12x12 [SIJD] 

18x12 [S**D] 

12xl8G[SD] 

6x6 

[S2D] 

12x12 [SID] 

18xl2G[SD] 

12x18 [VD] 



12x12 [S2D] 

18x12 [VD] 

12xl8G[V**D] 



12x12 [S*D] 

18xl2G[VD] 

12xl8G[VD] 



12x12 [S*ND] 
12x12 [VD] 

18xl2G[V**D] 




18x18 [SD] 

24x24 [SD] 




18x18 [SND] 
18xl8G[SD] 

18xl8G[SND] 

18x18 [VD] 

18x18 [VND] 
18xl8G[VD] 

24x24 [VD] 


where 





6x6 \ 

12x12 ) 

• f 





• 1 

12x18 ) 

m % 

number of rectangular subdomains 

in X-Y directions 


• 1 
• 1 

24x24 / 





G smoothly graded (as opposed to abruptly varying) 

S sinusoidal initial conditions with maximum u perturba¬ 

tion of 5.5 msec“^ and U=10 msec"^ 

S* same as S but one-fifth of the perturbation amplitude 

S** same as S but U=0 

V vortex initial conditions with maximum tangential wind, 

Vij, =25 msec“^ and U=10 

V** same as V but U=0 

1 diagonal slant from upper left to lower right 

2 diagonal slant from lower left to upper right 

(no number implies alternating slant) 

D diffusion included 

ND no diffusion 
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12 depending on the mesh. Longitudinally averaged wave amplitude infor¬ 
mation was obtained and used in data reduction. These amplitude data 
were used to measure generated noise and damping effects. Phase informa¬ 
tion was used to obtain wave speed by computing phase changes at two hour 
increments. 

The spectral analysis scheme required equally spaced data along 
a latitude circle. Since some grids tested had unevenly spaced nodal 
points, the interpolation property of FEM was used in these cases. If a 
value of a field variable, z, is known at an element's vertices located 
at points (x.,y^)(i=l,2,3); then the following is true: 

t 

= a + bx^ + cy^ 
z^ = a + bx^ + cy^ 
z^ = a + bx^ + cy^ 

Solving this 3x3 system for a,b, and c; one can obtain z at any 
point P interior to the element. 

3. Matrix Inversion Technique 

A straightforward, but time consuming, Gauss-Seidel technique was 
used for linear system solution for vectors of the three variables in 
each time step. It was noted that 10 to 15 passes were necessary for 
system convergence with no explicit boundary conditions. Eighteen to 23 
passes were required when boundary conditions were imposed. The condi¬ 
tions were explicitly imposed in the continuity and meridional momentum 
equations. Also of interest is the fact that convergence was only 
possible,when applying boundary conditions, if the direction of G-S 
iteration was reversed every pass. 
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Convergence was considered to be achieved when the following 


criterion was met: 


n n-ll 

X **• X I 

i i ‘ max over i 


< 1 X 10 


-6 



max over i 


with: n being a G-S iteration 


i being a nodal point number 


x^ being an element of the solution vector, {x} 


A more efficient scheme such as sequential over-relaxation is 


desirable to save time* Actual project completion required the use of 
many meshes all of which probably possess different optimum over-relaxa¬ 
tion factors, so G-S was used for expediency* 

4* Diffusion Coefficients 

A working maximum stable diffusive coefficient, which is 

equal to a comparable finite difference value was defined as: 



Using this as a maximum, various values of were tried on the 
12x12 grid. Values ranging from *0005 to *05 of k were tried and the 


spectral and visual results compared* The ideal was to use as little 
possible diffusion but enough to control generation of noise in higher 
wave numbers* A subjective value was obtained of *005 which was 
used in all tests including diffusion* The size of was limited by 
the minimum grid separation of each mesh* Values used are shown below: 
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TIME STEP MINIMUM GRID SPACING 


GRID 



(minutes) 

(km) 

/ 2 -1, 
(m sec ) 

6x6 



16 

580 

4.44x10^ 

12x12 



8 

290 

2.22x10^ 

18x12 

& 

18xl2G 



C 

12x18 

& 

12xl8G 

4 

145 

1.11x10 

18x18 

Sc 

18xl8G 




24x24 


4 




Final 

graded 

0.5 

8.26 

(see text) 


5. Data Reduction and Presentation 


The following noise or roughness measure was defined: 


Roughness = 


N 


i: 


a. 
1 


a 


1 


with a = zonally averaged amplitude from spectral 
analysis 

i = wave number index 

N = maximum Fourier wave number (3, 6, or 12) 


A measure of damping was used which was the ratio of amplitude 
of any Fourier component geopotential perturbation to its amplitude at 
t=0. This was used to determine the detrimental effects of diffusion 
on the solutions. 

An analytic wave speed was derived from an average propagation 
speed determined by phase angle changes. The quasi-geostrophic value 
of the phase speed, c, is given by 

(k + k - X ) 

- X _y_ 

^ ^ 2 2 2 

(k ^ + k ^ + X^) 

X y 
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a 


with 


k 

X 

k 

y 


= wave number in x direction 
= wave number in y direction 
= 


o 

(j)^ = mean free surface height 
U = 10 msec ^ 


V 


This formula gives phase speeds of 9.857, 9.958, 9.981, and 
9.989 meters per second for wave niombers 1, 2, 3, and 4 respectively. 

Graphical plots were produced for all tests. Those not dis¬ 
cussed as part of the results are presented in Figures 35 through 54 
inclusive. The results for the entire length of the cyclically con¬ 
nected channel were not plotted. Each plot extends to within one 
regular grid increment of the boundary where the mesh is connected. 

The maps of (p and u have the domain mean removed. The nodal points 
are marked with small crosses. 


B. SINUSOIDAL RESULTS 

The effects of diagonal slant can be seen in Figures 10 and 11. 

The most noticeable difference exists in Figure 10 where the mirror 
image effect of diagonal bias shows itself in the latitudinal direction. 
The cases shown possessing an alternating slant form a rough average 
between the opposite slant cases. Figure 11, which shows the (p ampli¬ 
tude, gives a slight hint of bias effects, but the tendencies are less 
pronounced. Wave speed data for similar cases on the 24x24 grid were 
not obtained, but one would expect less bias as the element size 
approaches zero. A plot of wave number one's amplitude versus latitude 
for the 12x12 mesh showed no significant variation for the three differ¬ 
ent diagonalization methods, so apparently the effect does diminish 
with finer resolution. 
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Figure 10. Latitudinal variation of wave speed (48 hour time averaged) 
for four grids. 
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Figure 11. Wave number one amplitude at t=48 hours. 






The latitudinal variation of phase propagation of various grids is 
shown in Figure 12 . In general all phase speeds are greater than 
analytic, but not significantly so. It should be noted that the scale 
of the abscissa in Figure 12 covers a very narrow range of phase speeds; 
even the coarsest grid is within four per cent of the true phase speed. > 
The phase velocities depart most from analytic near the boundaries, 
where the amplitude is small. A somewhat anomalous result is the 
better treatment of propagation by the 12x12 as opposed to the 24x24 
mesh. This is most apparent below the mid-channel latitude. No explana¬ 
tion is offered for this condition. 

The necessity of some spatial smoothing is apparent when one compares 
Figures 13 and 14; Figures 15 and 16; Figures 17 and 18; and Figures 19 
and 20. Even if it were not for this unwanted noise, large phase errors 
result without smoothing, as Figure 21 suggests. Even by decreasing the 
initial disturbance amplitude by a factor of five, the phase errors are 
almost identical with unsmoothed tests. 

The fine, course, and intermediate resolution grid results are shown 
for their roughness properties in Figure 22. Even with four times the 
damping, the 6x6 mesh cannot compete with the 24x24 mesh in this compari¬ 
son. The 6x6 results become even less desirable when one sees that the 
wave amplitude is cut more than in half (Figure 23). An enlarged view 
of Figure 22 is presented in Figure 24. This diagram should show if any 
significant improvement can be gained by using a grid which is inter- 
mediate between the 12x12 and 24x24 mesh, thus increasing accuracy with 
less computational sacrifices. It should be noted that all these inter¬ 
mediate grids possess one half the diffusion of the 12x12 mesh. 

Apparently two of the grids tried do offer the desired quality, the 12x18 
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Figure 12. Latitudinal variation of wave speed (48 hour 
time averaged). 
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Figure 13. 12x12 [SD] at t=48 hours. 
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Figure 14. 12x12 [SND] at t=48 hours. 
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Figure 15. 12x12 [S*D] at t=48 hours. 
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Figure 16. 12x12 [S*ND] at t=48 hours. 
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Figure 17. 18x18 [SD] at t=48 hours. 
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Figure 18. 18x18 [SND] at t=48 hours. 
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Figure 19. 18xl8G [SD] at t=48 hours. 
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Figure 20. 18xl8G [SND] at t=48 hours. 
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Figure 21. Latitudinal variation of wave speed (48 hour 
time averaged) with no diffusion. 
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Figure 22. Latitudinal variation of roughness for three grids at t=48 hours 




































































Mesh 

Figure 23. Damping characteristics of various grids at 
t=48 hours. 
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Roughness 

Pigure 24. Latitudinal variation of roughness at t=48 hours 

































and 12xl8G meshes. Both of these meshes have better resolution over 
the region where the disturbance is strongest. It was hoped that the 
grid with a graded rather than abrupt change in the y direction would 
be better. However, it appears the graded version (12xl8G) is only 
better at latitudes south of 15 North. The 18x12 and 18xl2G meshes 
are generally worse than the 12x12 coiinterpart although the 18xl2G is 
the better of the two. This perhaps is the result of the disturbance 
being advected across this graded or abruptly changing region of the 
two grids. One would, however, hope that the 18xl2G would show better 
noise properties than the 12x12. It is plausible that the element size 
change is still too abrupt to describe the scale of motion used with 
an advecting velocity. The remaining two cases to be identified are 
the 18x18 and 18xl8G. The 18x18 case shows an oscillating effect pro¬ 
ceeding northward. This could be explained as some pseudo-boundary 
influence induced by the abrupt element size change. The 18xl8G handles 
the situation somewhat more effectively. The oscillation is still 
present, slightly reduced, but the roughness is significantly less than 
the 18x18 case and is better than the 12x12 case at most latitudes. 

The spectral distribution of wave energy at 48 hours is shown in 

2 -2 

Figure 25. Some grids had disturbance amplitudes of less than 1 m sec 
(gpm) and therefore don't appear on this logarithmic plot. As expected, 
the 24x24 has the least tendency for energy transferral to higher wave 
numbers. Of interest is the large amount of energy exchanged by the 
18x12 grid. It has the largest energy levels in all higher wave numbers 
save wave number four, where its smoothly graded companion (12xl8G) is 
slightly higher. Relating the two 18x18 grids, one can see an expected 
improvement when a graded version of the same grid is used; a useful 
feature which is not as evident at latitudes different from 15. 
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Figure 25. Spectral characteristics with sinusoidal 
conditions at t=48 hours. 
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As an interesting sidelight, consider Figure 26. The gravity wave 
fluctuations due to slight imbalances are indicated on this diagram. 

The free surface gravity waves allowed in this case have an analytic 
period of 4.4 hours calculated for a wave length equal to that of the 
channel. This agrees well with the short period fluctuations shown in 
wave speed. 

C. VORTEX RESULTS 

The gravity wave influence on geopotential height variations is also 
implied in Figure 27. The diffusion is also evident as the stationary 
vortex fills through 48 hours. The gravity wave amplitude in this case 
is larger initially than with the sinusoidal case because of the larger 
tr;incation errors in representing the initial state. Rapid filling oc¬ 
curs in the first 18 hours due to diffusion. 

Again diffusion is necessary. To verify this compare Figures 28 
and 29. The undesirable requirement of significant diffusion at 500 mb 
is obvious, but in the vortex cases the diffusion imposed significantly 
danps the disturbance. Further experiments should verify if the scheme 
could survive with less diffusion of the type used. 

Figure 30 shows a view of damping for grids other than the 6x6 and 
24x24. Eighty eight per cent of the original mid-latitude perturbation 
is contained in waves one through three, so these should be of primary 
concern. All grids in this figure except the 12x12 have identical damp¬ 
ing coefficients and therefore should treat up to wave number three 
similarly. This is the case except both 12x18 grids don't damp wave 
nximber three as much. Since this Laplace type filter is not very wave 
number selective, one would expect higher diffusion of longer waves. 
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Figure 26. Temporal wave speed fluctuations. 
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Figure 27. Temporal geopotential fluctuations at center of stationary 
vortex. 
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Figure 28. 18x18 [VD] at t=48 hours. 
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Figure 29. 18x18 [VND] at t=48 hours. 
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Figiire 30.• Damping characteristics with vortex at t=48 hours. 
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It appears that both 12x18 meshes are not consistent with this assump¬ 
tion, Both damp wave number three about the same as two but treat wave 
numbers four and five in opposite senses. The 12xl8G grid also shows 
the same tendency as the 12xl8G though it damps wave number three more. 

The 18x12 and 18xl2G grids also show similar treatment of waves four 

✓ 

and five. The large damping of wave five relative to four can be ex¬ 
plained in part by the fact five's initial amplitude is about two and 
one half times that of four. The large damping of wave two by the 12x12 
mesh is a serious defect of the required amount of coarse grid diffu¬ 
sion. The actual spectral results of these grids are shown in Figure 
31. 

Figures 32 and 33 show the final distributions of amplitude for waves 
one and two, respectively. The large damping may make it hard to dis¬ 
tinguish between results for grids with greater resolution than the 
12x12. If this is but a secondary effect then one could say that extend¬ 
ing the entire domain resolution to that of the 24x24 is not necessary. 

A comparison of mid channel wave speeds shows all grids studied give 
very accurate results for wave numbers one and two. Phase propagation 
of wave number three is handled poorly by the 12x12 grid and is 270% 
in error for the 12x12 case and 70% in error for the 18x12 cases. All 
grids but the 24x24 are grossly in error for wave number four; at best 
175% error (18x18 grids) and at worst 240% error for the 18xl2G case. 

The 24x24 case is 25% in error for wave number four. All per cent errors 
mentioned are over estimates by the grids in wave speed. Phase predic¬ 
tion accuracy rapidly deteriorates for higher wave numbers. Propagation 
errors for large time integrations are contaminated by artificial noise 
being filtered to higher wave nimbers. 
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Figure 31. Spectral characteristics with vortex conditions at 
t=48 hours. 


73 























































































































74 


Figure 32. Latitudinal variation of geopotential for wave number one at 
t=48 hours. 
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Figure 33. Latitudinal variation of geopotential for wave number 
two at t=48 hours. 
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Studies with the highly variable grid shown in Figure 6 were cur¬ 
tailed. It was therefore not possible to fully investigate the movement 
of the grid with the vortex. A few experiments were carried out using 
both the vortex and single harmonic wave. The experiments with simpler 
initial conditions are not presented. A study was made of the vortex 
with no mean flow. The results of this investigation are shown in 
Figure 34. 

Basic tests were done with a constant over the domain. The value 
for this coefficient was again selected based on the smallest grid in¬ 
crement. On this grid, this value was very small, so small in fact 
that another formulation of diffusion was required. 

A more involved diffusion formulation was tried which assumed a 
diffusion coefficient constant over each element. A vector type diffu¬ 
sion was then applied. Simply stated, the diffusion over an element 
was derived using: 

■Se ' ^ 'Sm 

with: = diffusion on an element, e 

A = area of element e or element 1 

= diffusion coefficient based on minimum grid 
spacing, associated with elements 1 through 
24 (elements surrounding the nodal point in 
the center of diagram Figure 6). 

The plots in Figure 34 show that an initial roughness developed 
evident at six hours due to initial gradient imbalance which at later 
times was minimized as smoothing and geostrophic adjustment effects 
dominate. 


76 







K(lcm) 


<P Disturbance (gpm) 
6 hours 



0 Disturbance (gpm) 
24 hours 



0 Disturbance (gpm) 
48 hours 


Figure 34. Geopotential field for stationary vortex on final mesh. 
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This final case was integrated on both an IBM 360/67 machine and a 
CYBER 175 machine. The CYBER machine is twice as accurate so machine 
truncation effects were minimized. A comparison of spectral results 
frcm both machines at t=12 hours showed little difference, so machine 
truncation is not of concern. 
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Figure 36. 6x6 [SID] at t=48 hours. 
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Figure 37. 6x6 [S2D] at t=48 hours. 
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"Figure 38. 18x12 [S**D] at t=48 hours. 
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Figure 39. 12x12 [SID] at t=48 hours. 
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Figure 40. 12x12 [S2 d] at t=48 hours. 
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Figure 41. 18x12 [SD] at t=48 hours. 
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Figure 42. 18xl2G [SD] at t=48 hours. 
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Figure 43. 12x18 [SD] at t=48 hours. 
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Rigure 44. 12xl8G [SD] at t=48 hours. 
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Figure 45. 24x24 [SD] at t=48 hours. 
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Figure 46. 12x12 [VD] at t=48 hours. 
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t’igure 47. 18xl2G [V**D] at t=48 hours. 
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Figure 48. 18x12 [VD] at t=48 hours. 
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Figure 49. 18xl2G [VD] at t=48 hours. 
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12x18 [V**D] at t=48 hours. 
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Figure 51. 12x18 [VD] at t=48 hours. 
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Figure 52. 12xl8G [VD] at t=48 hours. 
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'Figure 53. 18xl8G [VD] at t=48 hours. 
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Figure 54, 24x24 [VD] at t=48 hours. 
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XI. CONCLUSIONS 


% 

In this research a finite element model has been developed for the 
barotropic primitive equations. The motion is confined in a periodic 
channel between rigid walls. The model was tested with two sets of analy¬ 
tic initial conditions. 

The experiments with a sinusoidal wave showed excellent phase pro¬ 
pagation when compared to the quasi-geostrophic analytic solution and 
when some diffusion was present. Even the 6x6 grid gave an error of only 
about 4 per cent which is no larger than the error obtained with a fourth 
order finite difference model. All grids produced computational noise 
without diffusion which made the forecasts unrealistic at t=48 hours. 

In order to control this noise Laplacian diffusion terms were added to 
the momentum equations. The low resolution experiments required a much 
higher diffusion coefficient than the high resolution experiments. This 
lead to heavy damping of the original wave for low resolution and little 
damping for high resolution. 

The same initial disturbance was tested with grids which contained 
an abrupt change in element size and grids which contained a gradual 
variation in element size. These variations in element size did generate 
some noise, but this noise was handled by the previously introduced 
diffusion. The smoothly varying elements gave somewhat better solutions. 

The model was also tested with an isolated vortex with the initial 
pressure determined analytically from the gradient wind equation. The 
experiments with a uniform element size displayed gravity waves which 
were excited by the finite element imbalance in the initial conditions. 
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The required diffusion dair5)ed and enlarged the initial vortex. The vortex 
experiments with a variable element size showed much better results. This 
was a result of the higher resolution around the vortex and the lower 
diffusion. Some experiments were performed with a continuously varying 
element formulation which had very small central elements. This procedure 
was very uneconomical because a very small time step had to be used 
throughout the domain (see Table II). 

This study has verified the previously discussed phase accuracy of 
the finite element method. However, the noise generation was larger than 
expected. It is not known whether or not some of the noise arises from 
imperfect boundary conditions at the walls. Cullen (1974, 1976) has 
discussed the necessity of high order smoothing with finite element models. 
In any case the model developed in this study would be benefited by a 
high-order spatial filter which would control the noise while not affect¬ 
ing the larger scale features. 

The main problem with the use of variable element size is the require¬ 
ment that the same time step be used everywhere. Finite difference 
models need use small time step only in the high resolution regions. 

Until this problem is solved it would appear that the element size should 
not be varied by a large amount. But since the FEM is more accurate 
perhaps tropical sto 2 nn forecasts can be made without excessively small 
element sizes. The model developed in this paper requires testing with 
a moving element system. Further work is also required to make the model 
more efficient. 
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TABLE II 


48 Hour PROGNOSIS 

IBM 360/67 TIME AND STORAGE REQUIREMENTS . 
(Fortran H compiler) 


GRID 


AVERAGE 

EXECUTION TIME 


STORAGE REQUIREMENTS 
(IN BYTES, 1 BYTE = 4 


6x6 

12x12 

3 minutes 

22 minutes 

49K 

129K 

12x18 

18x12 

1 hour 20 minutes 

188k 

18x18 

1 hour 50 minutes 

26 IK 

24x24 

3 hours 40 minutes 

44 2K 

Final 

Graded 

Grid 

10 hours (IBM) 

1 hour 20 minutes (CYBER 175) 

168K 


WORDS) 
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